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Abstract 



Ci • The Ewald summation technique is generalised to power-law l/jr"! potentials in three-, two- 

and one-dimensional geometries with explicit formulae for all the components of the sums. The 
^ ' cases of short-range, long-range and "marginal" interactions are treated separately. The jellium 



model, as a particular case of a charge-neutral system, is discussed and the explicit forms of the 
Ewald sums for such system are presented. A generalised form of the Ewald sums for a noncubic 



— Jj I (nonsquare) simulation cell for three- (two-) dimensional geometry is obtained and its possible field 

Q ' of application is discussed. A procedure for the optimisation of the involved parameters in actual 



simulations is developed and an example of its application is presented. 



I. INTRODUCTION 

The behaviour of many-body systems is often governed by the long-range Coulomb poten- 
tial between charged particles. Numerical simulations of such systems are usually performed 
by considering a finite number of particles in a cell with periodic boundary conditions. The 
correct estimation of the potential energy in such systems requires of a summation over all 
images created by the periodic boundary conditions. For long-range interaction potentials 
such direct summation either converges slowly or it is conditionally convergent, making its 
evaluation computationally cumbersome. Instead, the performance of the calculation can 
be greatly improved by using Ewald summation methods [l|. In these methods, the slowly 
convergent tail of the sum in the potential energy is represented by a rapidly convergent 
sum in momentum space. The method is named after Paul Peter Ewald who in his pio- 
neering work dated almost a century ago calculated the electrostatic energy in ionic crystals 
(a detailed derivation for the Ewald sums for the Coulomb potential can be found in the 



• B). 



work of de Leeuw et al. [2i). An alternative approach to deal with long-range systems is 
proposed by Smith [3|. In his method, the Hamiltonian and equations of motion are derived 
using constraints on the velocities of particles. Instead, in the following we will stick to a 
standard model for the Hamiltonian and will consider ways to improve the convergence in 
the potential energy. 

For a good performance in simulations of large A^-particle systems, a number of modi- 
fied summation methods has been developed. Historically, the first efforts to enhance the 
Ewald method consisted in looking for appropriate truncation schemes, but all of them were 
strongly dependent on the system properties, in particular on the system size. Tabulations 
of precalculated terms in both real space and momentum space sums [4| , as well as polyno- 



BJ? 



mial approximations of the involved functions [5|-[7|], were also proposed to look for a balance 
between calculation time and truncation errors. Nevertheless, these approximate methods 
suffer from error accumulation in simulations of large systems, and do not allow for reducing 
the overall 0{N'^) complexity of the original Ewald summation. The work of Perram et 
al. }p\ was the first to give a way to optimise the splitting of the interparticle potential be- 
twee,, the .o„g-range and short-range parts to yield a total complexity of O(N^n), a special 
modification of the Ewald method called Wolf summation [9|, [101], based on a damping of 
the Fourier-transformed part of the sum, was posteriorly developed in order to render the 



original Ewald summation more efficient for non-periodic systems and large model sizes. 

Another way for improving the Ewald method is to perform fast Fourier transform (FFT) 
of a reciprocal space sum on a mesh. The oldest algorithm of this kind is the so-called 
Particle-Particle Particle-Mesh (P^M) method, invented by Hockney and Eastwood in the 



late 80's ll|]. The P^M technique is based on a distribution of the charge density on a 
grid using a certain smooth assignment function and then the discrete Poisson equation is 
solved using FFT. This algorithm appeared to be less complex to yield 0{N In N) with an 
appropriate choice of the free parameters. The P^M algorithm was recently improved by Bal- 
lenegger et al. [14| for calculation of energies, bringing, as claimed, the maximal precision in 
the energy by an optimisation of the "influence" function (a substitution of the potential in 
the Fourier-transformed Poisson's equation). For a comprehensive introduction to Ewald- 
and mesh-based techniques we recommend to refer to the cited work of Ballenegger and 
coauthors where special attention is paid to the estimation of both sum truncation-imposed 



and grid- imposed errors. The extension of this method, called Particle Mesh Ewald [12 1 
(PME), makes use of the analytical form of the sum in the reciprocal space and evaluates 
potentials via FFT instead of interpolating them as P'^M does. Although PME is slightly 
more complex then the P'^M algorithm, it is still 0{N In A^) and allows to reduce significantly 
the memory expenses. Later Particle Mesh Ewald method was reformulated by Essmann 



et al. 13|, making use of cardinal 5-splines to interpolate structure factors. This approach, 
called Smooth Particle-Mesh Ewald (SPME) substantially improved the accuracy of PME 
with a comparable computational cost, as it still scales as 0{N\n.N). SPME is also claimed 
to be applicable to potentials of the polytropic form l/lrl'^. In general, the conventional 
FFT-based approaches suffer from the severe fallback of requiring equidistant particle posi- 
tions. The invention of the variant of Fourier transform for nonequispaced nodes (NFFT) 
opened a path to overcome this shortcoming, while keeping the introduced errors below the 
specified target levels. The nonequispaced fast Fourier transform is currently considered as 
a promising means to improve the Ewald summation performance, with open code imple- 
mentations availabl e | l5l| . The early variants of the NFFT algorithms are reviewed in the 
work of A. F. Ware Il6|: a general approach to the fast summation methods based on NFFT 

n 

can be found in the article of G. Steidl [17J. 

The most recent family of algorithms based on the Ewald approach are the tree-based 
algorithms, with the fast multipole method (FMM) being the most known and widely used 



among them. The algorithm, developed primarily by L. Greengard and V. Rokhlin [18|, is 
based on the idea of keeping the direct summation of potentials or forces for the nearby atoms 
and approximating the interactions of the distant atoms by their multipole expansions. FMM 
offers the asymptotically fastest performance among the Ewald-related algorithms, being 
linear in N in most cases and not worse than 0{N In N) with explicitly controlled accuracy. 
The FMM technique is naturally applicable to inhomogeneous and non-periodic systems, 
being also easy to parallelise since it is an entirely real-space summation. Since then the 
algorithm was significantly improved in efficiency, mostly by introducing new diagonal forms 
of translation operators [19]. However, FMM has an intrinsic shortcoming, when applied to 
molecular dynamics calculations, as the energy conservation it brings is poor; the method 
per se is also rather cumbersome in implementation. Another group of methods, based on 



the multigrid methods of solving elliptic (in this particular case - Poisson's) equations 



2Q|, 



was developed a decade ago [2l|. These methods allow to preserve the scaling 0{N) and 
parallelisation advantages of tree-based methods, as well as the applicability in simulations 
without PBC, being on the other hand satisfactorily energy-conserving and additionally 
accelerated on all length scales. 

A detailed comparison of the optimised 0{N^^'^) pure Ewald technique, FFT-based sum- 
mations, and multipole-based methods was made by H. G. Petersen |22| for systems with 
approximately uniform charge distributions, taking into account a possible parallel imple- 
mentation. According to Petersen, the method of choice with a number of particles below 
10^ is the conventional Ewald summation, PME is preferable in the range A^ ~ 10^ — 10^, 
and the fast multipole method should overperform them with A^ > 10^. A more recent and 
ample review of FMM, P^M and pure Ewald methods by Pollock and Glosli [23], based 
partially on their own calculations, implies that P^M is faster than the Ewald summation 
already for 500 particles, although it is stressed that the other factors as the ease of the 
coding, the system geometry, as well as the code optimisation can change the choice. We 
would also suggest a thorough survey of different Ewald summation techniques given in the 



work of Toukmaji and Board 24 1. 



An approach, alternative to using cubic periodic boundary conditions in a calculation of 
long-range interactions, called Isotropic Periodic Sum (IPS), was recently proposed by Wu 



and Brooks 



251 ]. The main goal of their approach is to deal with long-range interactions. 



avoiding artificial correlations and anisotropy bias induced by a PBG-based summation 



in a cubic box. In this technique, only the interactions of a particle A with the others 
within a certain radius Vc are taken into account (as in a plain cut-off scheme), and this 
spherical simulation zone is repeated in an infinite number of shifts by vectors rgh, such 
that |rsh| = 2A^|rc|. Therefore, the particle A interacts not only with B (within the sphere 
radius), but also with all the images of B, occupying homogeneously the shells of radii 
2N\rc\-, centered in B. The subsequent integration and summation over the shells allows 
to obtain explicit expressions of forces and energies for a number of interactions of most 
physical interest, like electrostatic, Lennard- Jones and exponential potentials. The method 
is known to yield a performance close to the one shown by the Ewald summation, but 
without imposing unwanted symmetry effects. 

Since its proposal, the Ewald method has been applied to a large number of physical 
problems, although mostly to systems with the Coulomb 1/|^| interaction potential. In a 
recent work by R. E. Johnson and S. Ranganathan 26|], a generalised approach to Ewald 
summation is stated to obtain potential energy and forces for systems with a power-law, 
Yukawa potential and electronic bilayer systems. The Ewald method for two-dimensional 
systems with electrostatic interactions was developed by Parry [27], but his technique ap- 
peared to be computationally inefficient. Spohr et al. [28!] studied a slab geometry by treating 
the simulation cell as a fully three-dimensional one with the conventional Ewald summation. 



Later on, a significant advance was made by Yeh and Berkowitz [29|, as the authors managed 
to obtain the explicit correction term for the rigorous three-dimensional Ewald summation, 
that brings the results for a slab system in a satisfactory agreement with the 2D summation. 



The 2D Ewald technique was also applied by Wen Yang et al. [30!] to calculate the energy of 
Coulomb particles in a slab system with a uniformly charged surface. Recent applications 
to dipolar bosons in a 2D geometry have been made by C. Mora et al. [31] and Xin Lu et 



al. 32| . On the other hand, the explicit forms of the Ewald sums for Yukawa interactions 



have been also reported: in 3D geometry, with partial periodical boundary conditions [33] 



and in 2D geometry 3j]. The Ewald method can also be useful even applied to fast de- 



caying power-law potentials. For instance. Shirts et al. in their recent work [35| argue the 
need for taking into account the effects of cutoffs in molecular dispersion interactions due 
to a Lennard- Jones potential, especially in non-isotropic and inhomogeneous media. The 
authors developed two formalisms for the estimation of these cutoff errors in binding free 
energy of macromolecular systems, which can in principle be extended to the other observ- 



ables. However, it is claimed that the adequate implementation of the Ewald summation 
for this kind of systems may render their corrections unnecessary by mostly eliminating the 
cutoff-dependent behaviour. 

In the present work, we report explicit expressions of the Ewald sums for the general 
case of particles interacting via a l/jrl'^ polytropic potential and in 3D, 2D, and ID ge- 
ometries. The closed derivation of these sums is given, with special attention being paid to 
conditionally convergent potentials. One of the difficulties of the derivation is that different 
terms have to be considered in the cases of short-range, long-range or "marginal" potentials. 
In the case of a short-range interaction, the original slowly convergent sum is represented 
as a linear combination of two rapidly convergent ones. For a long-range interaction, the 
condition of charge neutrality in the simulation cell is shown to be necessary to make the 
energy absolutely convergent within the considered scheme. The introduction of a uniform 
neutralizing charged background {jellium), as a particular case of a charge-neutral system, 
is also discussed. The explicit forms of the Ewald sums are reported for a jellium system and 
for an arbitrary polytropic potential. We explicitly calculate the expressions for physically 
relevant interactions as Coulomb, dipole-dipole, and Lennard- Jones potentials. Finally, we 
have extended the Ewald sums to the case of a noncubic simulation cell, that could be useful 
in simulations of hexagonal closed packed (hep) and two-dimensional triangular solids. In 
addition, the general derivation path given in this work may be used to obtain the forms of 
Ewald sums for other interaction potentials. 

The computational efficiency is another important issue of the practical implementation 
of the method. In fact, one needs to choose correctly a free parameter, appearing in the 
integral representation of the sums, and to decide which number of terms should be kept 
in spatial and momentum sums in order to reach the required accuracy. The choice of 
these three parameters affects the difference between the calculated result and the exact one 
as well as the calculation complexity. Therefore, a certain optimisation of the parameters 
is always required. In the present work, this optimisation process is formalised and it is 
shown that following the described procedure the overall computation time is significantly 
reduced. The accuracy of the result is shown to be kept under control, with the only cost 
of a preliminary benchmark calculation. 

The rest of the article is organised as follows. In Section [TTl we formulate the problem, 
develop the general Ewald approach and report explicit expressions for the Ewald sums for a 



polytropic potential in a three-dimensional cubic simulation cell. Sections IIIII and HVl contain 
derivations of the Ewald sums in two-dimensional and one-dimensional geometries, respec- 
tively. In Section |Vl the case of a simulation cell with different side lengths is considered 
for three- and two-dimensional systems. The final general expressions and their particu- 
larization to the most physically relevant cases are presented in Section I VII The practical 
algorithm for the parameter optimisation and an actual application of the Ewald method is 
discussed in Section IVIII Summary and conclusions are drawn in Section IVIIII 

II. EWALD SUM FOR AN ARBITRARY POLYTROPIC POTENTIAL l/\r\^ IN 
3D GEOMETRY 

A. Basic assumptions and initial sums 

We consider a system of A^ particles inside a cubic simulation cell of size L with periodic 
boundary conditions. Thus, each particle with coordinates r in the initial cell has an infinite 
number of images r + nL in the adjacent cells. The total potential energy is estimated by 

N N 

(1) 



n^E' 



2 






where 0(r) is the interparticle potential, rjj = r^ — Tj, and the prime in the first sum means 
that the summation over an integer vector n must be done omitting the term n = when 
i =j. 

B. Analytic development 

In many physical situations, the interaction potential between two particles i and j has 
the power-law form qiqj/\r\'^ with positive k and g^, qj being the generalised charges of the 
particles. This sort of interaction is generally referred to as polytropic potential. 

First, let us consider the case of short-range potentials. A; < 3. As we will see later, the 
potentials corresponding to fc > 3 give a similar result. For A; < 3, the right-hand part of 
Eq. ([T]) diverges and it can be made convergent only if the restriction of charge neutrality 



is required, i.e., when X]i=i ?« = 0- ^^ has also been shown [36| that for a pure electrostatic 



interaction the total energy ([T]) can be conditionally convergent even in a neutral simulation 



cell because of a higher multipole contribution. The energy and forces are therefore depen- 
dent on the order of summation, which can also be implicitly set by a choice of a convergence 
factor. The ambiguity usually appears in a form of a constant or a position-dependent term, 
vanishing in the limit L — )■ oo. Hence, the preference in one or another factor should be dic- 
tated either by physical properties of a particular system or by arguments regarding rates 
of convergence to the thermodynamic limit. For a general discussion on the convergence 



issues appearing in periodic boundary conditions, see Ref . [37| • The main idea of the Ewald 



g 



is to 



summation technique in the approach proposed by de Leeuw, Perram, and Smith 

2 

multiply each component of the sum by the dimensionless factor e"**" , with s > being 
a dimensionless regularizing parameter, making the sum absolutely convergent. Then, the 
limit s — 7- is taken, so that the singularity in the initial sum ([T]) can be explicitly separated 
into a term depending only on s, that finally can be cancelled due to the charge neutral- 
ity condition. We take a similar multiplier c{n, r, s) = e~*l"'+''l yielding the same rate of 
convergence (since < r < 1 in units of L). As the sum, multiplied by c, is invariant to 
an arbitrary substitution r — t- n + r, the chosen convergence factor allows to preserve the 
periodicity of the potential in order to avoid any possible artefacts in the final results. 

For the sake of clearness of the derivation, it is convenient to use reduced length units, 
that is to use the size of the box L as unity of length and substitute r^j by rijL. From now 
on, and for simplicity, we use the notation r for Vij and, in case of possible ambiguity, we 
will stick to the standard notation Vij. Also, we rewrite the potential energy by splitting 
the total sum ([1]) into two terms: Jqi (the sum of the interactions between a particle with 
all the other particles in the box), and /qo (the sum of the interaction of a particle with its 
own images, comprised of the components i = j in Eq. ([T])). Explicitly, 

n = -l(/oi + /oo), (2) 

with 



/o.=5: 



7,3 



E 



.l<i<j<N ' 'J ' 



n£: 



1 V- e— ^ ^ 



(3) 



^00 = ^ E V5:^^ (4) 



where the shorthand notation n = \n\ is used. 



First, let us focus on the /qi term, which we rewrite as 

l<i<j<N 

where we have defined the "screened" interaction potential "ipir, s) = ^^ g-«k+"'l /\r + n|^, 
extended from a single cell to the whole coordinate space. Since the total potential energy 
consists of a sum of pair interaction components, we may consider a single pair without any 
loss of generality. 

Let us apply the equation 

x-^' = —— / f-^e-'^' dt , (6) 

r(s) Jo 

representing the definition of the gamma-function, to the polytropic potential \r + n|~^. 
Then the function tp may be represented in an integral form, 

1 /""^ 

^(r, .) = ^^ y^ tt-1 Yl e-*l^+-l^e-^l^+-l^ dt . (7) 

We expect that the integral i^^ contains a singularity that will be located in the vicinity of 
zero. Therefore, we split this integral into two domains [0, a^] and [a^, oo), the corresponding 
integrals being denoted as ipf^^ and ipmi, where a is some arbitrary positive constant, 

ipir, s) = ifJunir, s) + ^inf (r, s) . (8) 

In the following, we analyze the two terms of the previous sum ([8]). 

1. The explicit analytical form of the term ipini{r, s) can be found 

r(|)^A2 ^\r + n\'^ r(|) 

(9) 
where T{a, z) is the incomplete gamma function. From the large distance asymptotic 
expansion of this function, one obtains that the above lattice sum is absolutely and 
uniformly convergent if s > and a > 0. Therefore, one may simply take the limit of 
vanishing screening s — ;■ 0, 

, . .s^ 1 ^ r(|,aV + np) 

^ 2' n ' ' 



2. The calculation of '?/'fin(^, s) is done by making a separate analysis of the n = case, 



,n^0 



^fin(r,.)=Cr(^,^)+C„="(r',3). 



(11) 



Explicitly, 



,717^0 



^fin (^>S) 



V'fiT^r-,.) 



3 

7r2 



^2^ ri^O 



t2 



-1 






3-exp 

2 



-n^n^ 



t + s 



+ 2'Kinr 



dt 



3 

7r2 



.2 il-i 



rft)7o (t + s 



dt , 



(12) 
(13) 



where we have used the Jacobi transformation 

3/2 



38, 


39 



y^g-s|n+r|2 



( - j N exp[-7r^n^/s + 2ninr] for n e Z^ 



applied to 



exp[— s|n + rp — t\n + rp] = exp[— (s + t)\n + rp]. 



(14) 



(15) 



We evaluate the integral ip'^^ (r, s) by the following analysis. Consider separately the 
following factor of the integrated expression from flT2|) 



M 



exp 



t+s 



{t + s 



(16) 



It is clearly continuous and bounded on (0, +oo) as a function of (t+s), also notice that 
^fc/2-i ig absolutely integrable on (0, a^) for k > 0. In accordance with the standard 
convergence test for improper integrals, the integral tp^^ (r, s) converges absolutely 
and uniformly with s being considered as a parameter. Then, the limit s — )■ may be 
carried out and the integral becomes 



,n^0 



i^Z^'ir^s) 



3 

7r2 



E 



,27riTir 









7r2 cos(27rnrl 



r(|) 



fc-5 

t 2 exp 



2 



TT^n^ 



TT^n^ 



dt 



a^ 



(17) 
(18) 



The function En{z) is the exponential integral function, and we have cancelled the 
imaginary part of the sum (IT7|) by grouping the pairs with n and — n. 
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Now, we analyze the second term of ip&nif, s] 



^r'ir,s) 



3 



i2 ,i;_i 

t2 



r(f)io {t + s\ 



dt . 



(19) 



In terms of a new variable v = s/{t + s), 



t 



n=0( 

fin I 



3 

772 



:i-v)2 



-1 



,(fc-3)/2^^^ 



(20) 



s/(a'^+s) V 2 



The integration of i/j'^^^{r, s) for a l/|r|'^ interaction has to be carefully analyzed as 
a function of fc: 1 < k < 3, long-range potential; k = 3, marginal case; and k > 3, 
short-range potential. 



(a) Suppose 1 < A; < 3. The resulting integral, 



^sr"(r',s) = Cn=>,^) 



3 

7r2 

rig) Js/{a'^+s) 



V 2 



fc-i 
V 2 



,(k-3)/2^^ (21) 



may be given explicitly in terms of incomplete beta- and incomplete gamma- 
functions. Expanding the resulting function for small s, 



V^, 



n=Oi 

fin > 



r,s 



k-3 



3-k 



2TC2a 



k-3 



{k - 3)r [|] 



+ 0{s) 



(22) 



It is easily seen, that the only divergent term in the expansion (122|) is the first 
one, which we define as 

'3-k' 



Sis) 



k-3 

s^2ttT 



(23) 



We remind that the choice of a convergence factor (that explicitly affects the 
summation order) may in principle lead to additional contributions in the total 
energy if the convergence of the sum is conditional (like for a charge-neutral cell of 
Coulomb particles with non-zero total dipole moment). In the original derivation 
of de Leeuw et al. [2|, the factor exp(— sn^) results in an additional dipole-like 
component in V'firT^' which breaks the periodicity of the potential and therefore 
complicates its use in simulations with periodic boundary conditions. Moreover, 
this procedure |2| yields a nonvanishing dipole term exclusively for fc = 1 in 3D 
geometry, with the rest of the sums remaining unchanged. From our point of 
view, this discontinuity points out to an nonphysical character of the dipole term 



11 



appearing in the case of the Coulomb potentiaL Nevertheless, in a number of 



studies 



36 



37l | it is considered as a first order correction when the convergence 



to the thermodynamic limit is analyzed. The mere fact that the results for the 
two different convergence multipliers coincide when fc > 1 is a consequence of the 
absolute convergence of the higher multipole contributions in this case. 

(b) Suppose k = 3. In this marginal case, the expression (12T|) may be integrated 
directly to yield the following logarithmic dependence 

7r2 / —2as — 2a^ 



^': 



n=0/ 

fin 



r, s 



r(f) 



(a^ + 5)2 



+ ln(s + 2^2 + 2aVs + a2) - In s 



that close to s = expands as 



Cir°(r, s) = -27T In s - 47r + 47r ln(2a) + 0{s In s) 



(24) 



(25) 



with the diverging term 

S{s) = -2TT\n s . (26) 

k 

(c) Consider the remaining option A; > 3. In this case, (1 — f )2~ is bounded from 
above and {k — l)/2 > 1. It means that the integral converges absolutely and 
the only finite contribution to the integral comes from the first (constant) term 
of the integral expansion for small s, 



Cn=°( 



r,s) 



2T[2a 



k-3 



(27) 



{k - 3)r [|] 

The second term of the total potential energy, Jqo ([2]) can be derived in a similar form to 
the first one. The procedure to find the form of '?/'(r', s) is repeated here with Vij = 0, hence 
the results are obtained straightforwardly via ( ITOl) . ( fTSj) . ( l25l) and ( !27|l . 

1 



N 



'00 



-T."' 



i=l 



^2,' n 



^T{l,a^n') 



w^ 



+ E 



n^O ^ 2 



T—O. JZ/ k-1 



TT^ra^ 



a^ 



a 



nl + r 



+ Cn=V,^) 



(28) 



with the term ip'^^ °(r, s) depending on the potential parameter k via ( 122|) . (125|) or ( 127|) . 
Putting all together, the potential energy can be written in a more compact form as. 



n 



= y^(-^oi + -^00 ) 



1 1^1 1 

i<j 1=1 i<j i 



(29) 
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with the generahsed potential, 






(30) 



A constant shift in the definition of ip is introduced to satisfy by the property /^^^ ipdr = 0, 
convenient for a proper treatment of the background contributions (see Appendix). The 
functions entering in Eq. (!30|) are defined as 



R{n,r) 



r(|,aV + nP) 



K{n, r) = K(n) cos{27inr) 



(31) 
(32) 
(33) 



with 



Kin) = ; Ek-1 



Tr'^n'^ 



a^ 



The exphcit form of the function S{s) depends on the k value, 



S{s) 






s 2 2nT [\^] if fc < 3 (singular term) 
— 27r In s if A; = 3 (singular term) 

if A; > 3 



and the term ^ depends only on the choice of a, 



^ = J2{p{n) + K{n)) + C, + C2 



with 



p{n) 



r(|,aV) 



and K,{n) defined in Eq. ( 1M|) . The constants Ci and C2 are explicitly. 



(34) 



(35) 



(36) 



(37) 



An + An ln(2a) if A; = 3 
C2 



a 



r(f + i: 



(38) 

(39) 
(40) 
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C. Removing singularities for k < 3 

The diverging part lis (containing a singularity) of the total potential energy equals to 

n. = ^ E ^^^.■^(^) + ^ E ^'^(^) = ^ f E ^) ^(^) (41) 

i<j i \ i / 

and vanishes, if the charge neutrality condition ^^ gj = is taken. 

Consider now a charge-neutral system with a neutralizing background consisting of a large 
number of identical uniformly distributed particles of the opposite charge (the "jellium" 
model). We denote the numbers of negatively charged particles g_ and positively charged 
(background) particles g+ as N_ and A^+, respectively. By imposing charge neutrality, 
g+ = —[N_/N^]q_, with N the total number of particles, A^ = A^_ + A^+. 

The potential energy for the jellium model can be written as 

n = ;^ E mM-r^ol^) + %^. H ■ (42) 

The second term in Eq. (H2!) has a component proportional to N^q]^. Note that the 
negative charges g_ and their number A^_ is defined by the problem and therefore fixed. 
Hence, in the limit A^+ — t- oo, this term cancels Nj^q^ = {N^q'^)/N^ — )■ 0, and therefore this 
background contribution may be eliminated to yield 

Nn"^ + N+ql N_q^ , , 

Concerning the first term of Eq. (H2|) . let us split it into three pieces, 

\ E 9.g,^(^) = ^(5- + 25_+ + 5++), (44) 

\<i<j<N 

where the first sum corresponds to the interaction between the negative charges 

s- = E ^^^^•^(^) ' (45) 

l<i<j<N- 

the second sum is the interaction of the negatively charged particles with the positive charges 
of the background 

N- N++N- 

5-+ = E E m^i^)^ (46) 

4=1 j=l+7V_ 
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and the third one is the interaction between the background charges 

^++= E mMr)- (47) 

l+N-<i<j<N++N- 

The last two terms S ^ and 5++ are easily shown to be zero in the limit A^+ — )■ oo as a 

consequence of the zero value of the integral of ip over the simulation cell (see Appendix). 

With the above considerations we can finally write the expression for the potential energy 
within the jellium model as 

^''' = %Y.^(r,^lL) + ^i (48) 

In the more general case of different charges in a charge-neutral simulation cell (with 

a long-range potential) or a system with an arbitrary short-range potential the potential 

energy is given by 

1 sr^N J2 

n^^^' = Z^ E <i^<iMn,/L) + M^i ■ (49) 

A certain analytical conversion of the sum in the reciprocal space is also possible in order 
to sum it up faster. Expanding the sum that defines K{n^ r) (1521) . one can simplify it in the 
following way, 

X] ^i^i 5Z ^*^^' "^^ = 5Z '^(^) 5Z ^*^J cos(27rnr) 

i<j n^O n^O i<j 



1 . I 

- 2^ K(n) 2^ Qigj [cos(27rnrj) cos(27rnrj) + sin(27mrj) sin(27rnrj)] N, Q^ /, 't(n) 

2 

-^K{n) ^qjexp{2mnrj) - ^ J^g^^ J^ /«(n) (50) 



2 , 



J 



In this form, the sum over all pairs of particles in the reciprocal space is represented as 
a single sum over particles and thus it scales as 0{N) instead of (9(A^^). Notice that the 
number of prefactors K,{n) and exponents in the sum depends on a chosen cutoff, which in 
general also might depend on A^, making the overall complexity of the fc-space grow. Naive 
schemes with a and the cutoff not depending on A^ do not take into account the interplay 
between the r-space and fc-space sum complexities, thus leaving at least 0{N'^) in one of 
them. Nevertheless, as we show later, optimisation with a and cutoff depending on A^ 
gives a best total complexity of 0{N^^'^). An alternative method to sum up the momentum 
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space part is to use Fast Fourier transform-based techniques (like PME), which is fast as 
O(A^lnA^). 

The last term in Eq. fl50l) cancels the K{n) component of ^. Introduce the notation, 

^(r-)=^i?(n,r) + Ci (51) 

n 

i=J2pin) + C, + C2 (52) 

5'equai(^)= Q- "^ exp {2ninr j / L) (53) 

j 
Sq{n)= y^ Qj exp{2mnrj/L) , (54) 

j 

where S'equai is used when the system of equally charged particles g_ is considered. Within 
this notation the potential energy may be rewritten in the following forms, which are more 
efficient for numerical implementation, 

n^'' = § E ^(^^^V^) + ^ E <n)\S^.Un)\' + ^i (55) 

i<j n^O 



■Kj 

with Tj, Tij in the original length units. 



D. Short-range potentials and the marginal case 

In case of a short-range interaction (A; > 3), the potential energy does not diverge, which 
is clear from the form of the singular term S'(s)(see Eq. [35|l . Hence, there is no need to add 
a neutralizing background and, even more, the background must be necessarily excluded 
since it leads to a divergence in the energy. This is easily seen by considering the potential 
energy of the background separately 

/■cell ^ 

n.. = ^/ Uh^' (57) 

that contains a singularity in zero. The expression for the potential energy is simply equal 
to Eq. (iHD, 

n = 71 E mMr.,m + H^i . (58) 
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When k = 3 (marginal case), both ultraviolet and infrared divergences arise in zero for 
the background as well as in the vicinity of infinity (the logarithmic divergence in the energy 
of negative charges). The only coherent model here is a plain "quasi-neutral" gas consisting 
of a mixture of a finite number of charges per box with the constraint "^Qi = 0, i.e., with 
the positive background excluded. 

III. EWALD METHOD FOR TWO-DIMENSIONAL SYSTEMS 
A. General notes for lower dimensions 

The Ewald sums can be extended to two-dimensional (2D) systems interacting through 
polytropic potentials. The difference with the 3D case comes from a different form of the 
Jacobi imaginary transformation for the Jacobi 6'-functions [its 3D form is given in Eq. (JHj)]. 

The "third" Jacobi 6'-function Os^z^t) is defined as 

+00 
^3(^|r) = Y, e^"™'e2™" , (59) 

n=— 00 

and satisfies the Jacobi imaginary transformation, 

e,{z\T) = {-\rr^l^^^"'^'l- QAzT'\r') , (60) 

with r' = — 1/t. Under the change of variables, z = vrr and r = ivr/s, the ^-function 
becomes a Gaussian, which is the relevant function for performing the Ewald sums, 

+00 +00 

n=— 00 n=— 00 

This expression will be used later, in the derivation of the Ewald sum in one-dimensional 
systems. Equation (1^ may be easily generalised to the 2D geometry, 

^ e-'l^+^l' = (tt/s) ^ e--^"V^e2™r . (62) 

n n 

Comparing this result for 2D with its ID (!6T]) and 3D (TT^ counterparts one finds that the 
dimensionality D affects only the constant multiplier as (tt/s)^/^. 
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B. Derivation 

The analytical derivation of the Ewald sum in 2D proceeds similarly to the one already 
presented for 3D. Equations from ([2]) to fITT]) are also valid here because their derivation 
is done without explicit reference to the dimensionality of the problem. In particular, the 
integral ipini{r, s) converges absolutely and to the same value 



'4'ini(,r,s) 



s^O 



—y 



r(|,a2|r + n\ 



n\ 



(63) 



We make the same decomposition of the integral ^fin(^, s) as in 3D, 



with 



^L^'ir,s) 



71 



E 



^(2^n^0"^0 



t2- 



it + s] 



exp 



— TT^n^ 



t + s 



+ 27rinr 



dt 



^, 



n=0, 

fin > 



TT 



r,s 



U^ 



r(f)7o it + s\ 



■ dt 



(64) 



(65) 



(66) 



where the two-dimensional variant of the Jacobi transformation fl62l) is used. The difference 
between the pair of equations (!65| [66]) and their three-dimensional analogues ( !T2t [T3|) relies 
in a substitution of the 3D factor (7r/(t + s))'^/^ by the 2D one 7r/(t + s). 

First, we consider the term ■j/'gjf (r, s). Following the same analysis as for its 3D coun- 
terpart, it can be shown that this parametric integral also converges absolutely. It yields 



t 



fin ' 



TT 



r,s) 



E' 



,27rin.r 






t2 exp 



vr^ri^ 



dt 



-^7ra)s(27rnr) fc_2p 
Z_^ r(k\ 2 



Ti^O ^2' 

/,n=0 



TT^n^ 



a^ 



(67) 



The modification of the integral "i/'firT i^ ^^^^ straightforward, since it requires specific 
integrations and expansions in the series for small s. Namely, we have to evaluate the 
integral 



J,n=0 
rfin 



IT 



r(|) 



U 2 



s/(a'^+s) 



k 
f 2 



-s'^/'-'dv 



(68) 



which is the 2D equivalent of Eq. (I2UI) . 

In the following, we consider separately the cases of long-range potential (1 < A; < 2), 
marginal interaction [k = 2) and short-range potential (A; > 2). 



18 



1. 1 < k < 2. As in 3D, the integral can be found analytically via the incomplete beta- 
and incomplete gamma-function with known series expansions for small s. Omitting 
these unnecessary intermediate expressions, we give the final expansion for ip^^^, 

z^^=o = s"-^ ,/' ,,, + , ^""""r' 1 + 0{s'/') . (69) 

^^" sin(f)r(|) (A;-2)r[|] ^ ' ^ ^ 



The first term of the expansion. 



2 

fc-2 7T 



S{s) = S— 7^^ TTT , (70) 

sin(f)r(|) ^ ^ 

clearly diverges when s — t- 0. Similarly to the 3D case, this term is cancelled in a 

charge-neutral cell and hence, 

"^''^ "(A:-2)r[|] ^^^^ 

2. k = 2. The integration of Eq. (l68l) is performed to yield in the limit s — ;■ a marginal 
logarithmic dependence, 

^^=0 = -TT In s + 27r In a + 0{s In s) . (72) 

As for the 3D geometry, the jellium model is inapplicable in this particular case since 
the energy of the continuous background diverges. Nonetheless the diverging compo- 
nent 

^(s) = -7rlns (73) 

can be removed if we consider a charge-neutral system with a finite number of charges. 
In this case, 

^^=o = 27rlna. (74) 



3. k > 2. The integral (1681) can be evaluated by taking s = 0, since its convergence is 
absolute, 

^'" {k-2)r[!l\ ^^^^ 

The second potential energy component, /qo ([2]), is calculated as in the 3D case. The 
result for 2D is 

N 



loois) = J2 <lKi^^t\0, s) + ^i„f(0, s) - Cr°(0, s)) (76) 

N 



j=i 



.r r(|K ^0^(1) H «^ J r(| + i)+^«'^ 
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C. Final expressions 

With respect to the 3D case, the changes in the 2D Ewald sum appear in those terms in 
which the Jacobi transformation is used, that is in «:(n) and Ci, 

Ci = V^r° (78) 

The other terms, namely R{r, n), p{n) and C2, are not affected by dimensionahty and may 
be taken directly from the previous section. 

Within the jellium model for a long-range potential {k < 2), the Ewald sum is given by 

n-' = §E^(-^./^) + ^^ (79) 

i<j 

A more general form, applicable to any system with a short-range potential (A; > 2), a 
charge- neutral system with long-range interaction (fc < 2), or a marginal {k = 2) potential 
is expressed as 

1 f ^ 

^"^^ = z^ E m^^^j/L) + ^ E ^' • (80) 

In the same way as for the 3D systems we can modify the sum in the reciprocal space, 
and with the same notations ( ISTl) - (15^ (p, R and the constants Ci, C2 are the new ones, 
corresponding to 2D case) the potential energy may be given by 

n^^= §E^(^^./^) + ^E'^(^)i^-uaip + ^e (81) 

n''^"^= Z^E9^9^-^(^'^/^) + ^E'^(^)l^^l' + ^^~- (82) 



IV. EWALD METHOD FOR ONE-DIMENSIONAL SYSTEMS 

As it has been commented before for the 2D case, the differences due to dimensionality 
are caused by the form of the Jacobi imaginary transformation. In the derivation for ID, 
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one needs the following ones 



1 /""^ 

- ' '^^ ' (83) 



(84) 



X 


; t e at 

r(s) Jo 


+ 00 


+00 


n=— oo 


n=—oo 


■oo 


+00 



n=— oo ra=— oo 

Similarly to what discussed in the previous section, the only terms to be changed are those 
where the Jacobi transformation is used, namely if)^ (in Jqi in a radial- dependent form, in 
/oo for r = 0). The difference arises from a different power exponent (1/2) in (I84p and (ISSjl . 
that is in (TT8|) A; has to be substituted by fc + 2 (and vr^/^ - by vr^/^, respectively), yielding 

As far as the term ip^^^ is concerned, we should perform a simple integration and do a 
series expansion for small s. 



«r°^^/' ^i#d. (87) 

The estimation of this integral depends on the k value. In the following, we detail this 
analysis. 

1. k = 1, the marginal case, 

^fiT° = ^(-lns-2 + 21n(2a)) + 0(.). (88) 

As before, we keep only the constant term, considering the diverging term absent due 
to the charge neutrality condition. Therefore, with r(l/2) = ^/^^^ one has 

V^^r° = -2 + 21n(2«) (89) 

2. A; > 1, the short-range potential, 

1/2 O^fc-1 

K:' = f(iy 13T + ^(^) + 0{s^'-'^^' In .) . (90) 

In the limit s — )■ 0, it yields 

resembling the 3D result (!27|) . with the change k ^ k + 2 (except in the F term) and 
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The final result for the one-dimensional Ewald summation reads 

^(r) = ^i?(n,r) + ^7r(n,r)) + Ci (92) 

e = $^(pH + «(n)) + Ci + C2, (93) 

where Ci = 4'J~^ is taken from the expressions (189|) (if A; = 1) or (19T|) (if A; > 1). 

For k = 1, the only consistent system is the charge- neutral one with a finite number of 
particles. In this case and for a short-range potential {k > 1) one the potential energy is 
given by 

n^en ^_J2 q,q^^{r,^/L) + ^^^ • (94) 

i<j 

Although the Ewald method is applicable to one-dimensional problems, there is a direct 
way to calculate the sums for polytropic potentials 

-. n=+oo ^ 

J^K ^-^ U- _|_ ^ ft 

n=— oo ' ' 

For k > 1, this sum can be represented as a linear combination of the Hurwitz zeta functions, 

71 E ^-^ = 7l(^^(^) + ^^(1 - ^)) • (96) 

_^ft Z—^ U- _|_ ij, ft _^ft 

n=— oo ' ' 

In particular, for /c = 2 the sum converts into a familiar expression used in the Calogero- 



Sutherland model 40|, |4l| , 



1 +°^ 1 2 

1 >r^ 1 TT^ 

Z2 Z^ |r + n|2 " L2sin2(7rr) ' (^ ^ 

n=— oo ' ' ^ ' 

Notice that the sum ( l96ll may be expressed in terms of trigonometric functions only for even 
values of k via {k — 2) times differentiation of Eq. ( P7|) . Anyway, the possibility to find exact 
expressions for infinite sums in ID suggests that the use of the Ewald method might not be 
needed, but we keep it as a possibly useful mathematical relation and for completeness. 

V. EWALD METHOD IN A RECTANGULAR BOX OF ARBITRARY SIDE 
LENGTHS 

A. 3D case 

A special and interesting situation arises if we consider a simulation cell in a more general 
way, as a rectangular box with different side lengths (L^., Ly, L^ in the corresponding di- 
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mensions). The need to deal with a box of unequal size lengths may occur in the simulation 
of a solid with a noncubic lattice (the simplest examples include a hexagonal closed packed 
crystal in 3D geometry), since the lattice vectors n in the sum over images on ([T]) are no 
longer orthogonal. Focusing our analysis to a 3D geometry, the potential energy is now 
given by 

1 r Af Af 

n = o E ' E E <^(^'^- + ^0^^) , (98) 



nci.& 



i=l j=l 



with fir = {rixLx + nyLy + UzLz)/ Lq, n^.,^,^ being integer vectors along the corresponding 
axis x, y, z. We have introduced the geometric average Lq = {L^LyLzY/^ and we will 
use reduced Lq units for rjj, and hence r^j will be adimensional. Repeating the standard 
procedure, we multiply the potential energy by a Gaussian term exp(— s|nr + rp) and, at 
the end, we take the limit s — > 0, separating the converging part, if present. We group 
separately the interaction with images of other particles Jqi and the interaction of a particle 
with its own images /qo, 



n = -^(/oi + /oo) 



(99) 



where 



'01 



'00 



E 



E 



I I \k 

n&I? \.l<i<j<N ' '■? '^' 



E xr^E^' 



2 ■' — ' in.,. 



(100) 
(101) 



i=l 



Comparing the relations ( 199|) - (llOip to the cubic case (E]) - (jl]), one notices that these 
relations remain unchanged if n is formally substituted by n^, and the constant coefficient 
1/L^ is replaced by I/Lq- Therefore, all the results found without the Jacobi transformation 
(IT^ remain the same with n^ instead of n. In particular, Eq. flTOj) transforms into the 
following 



^ 



inf 



m: 



E 



n.,. 



(102) 
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The Jacobi transformation (H^ in a noncubic box has the following form 



E' 



s|Tir+r|^ 



\ ^ ^-s{niLi/Lo+riy 






i=x,y,z rii 






n ( ' V'^ 

_'i=x,y,z ^ ^ ' "' ' 


i=x,y,z rii ^ 


vr^n? 


s{L,/Loy 



exp{27imiriLo/Li 



(tt/s)^/^ Y^exp(-7r^|nfcp/s) exp(27rinfcr) 



(103) 



"fe 



with Uk = n^Lo/L^ + UyLo/Ly + UzLq/Lz the normalised displacement vector in momen- 
tum space. The last equation is obtained from the original expression (TT4|) by a formal 
substitution of the vector n by n^. 

In order to calculate ipf^^ we first modify Eq. fllSp . 



exp[— s|nr + r|^ — t|nr + r|^] = exp [— (s + t)\nr + r|^] 



(104) 



then insert it into the relation (I103p . and finally separate the summand n = 0, 



^: 



3 ra^ 4^- 



fin 



ni)±:Jo {t+s)i 



r^t' + Cn=° 



exp 



t + s 



+ 2ninkr 



dt + 



TTZ 



t2- 



r(f)yo {t + s)l 



■dt 



(105) 



The subsequent derivation follows exactly the derivation for a cubic box, with the change of 
n by rij. and Uk for sums in the real and momentum spaces, respectively. The final result 
for a 3D system in a noncubic box can be summarised as follows 



^y- r(fc/2,a>, + rp) ,^ 



r(A;/2)|n^ + r| 



^^ r(A;/2,a>,p) 



vrf«-3,,,(2.^^,)^^^^ ^^rW^, ^C, (106) 



rifcT^O 



r(fc/2) 



«^ 



J2 



r(A;/2)|n^|'= 






,,fc-3 



™^ -i^..f^^i^~i+C,+C. 



r(fc/2) 



a" 



n=|iE^(^^^/^o) + ^^ 



(107) 
(108) 



«<J 



with the constants C\ and C2 defined in (138!) and (l39l) . As it was done in the cubic box, the 
potential energy may also be given with the momentum space sum (linear in A^). Applying 



24 



the definitions, similar to Eqs (15T]) - ([5i 

7/'(r)=^i?K,r) + Ci (109) 

rir 

i= 5^pK) + Ci + C2 (110) 

'S'cquai(nfc)= g- ^ exp(27rinfcrj/L) (111) 

i 
Sq{nk)= y^^qjexp{27iinkrj/L) , (112) 

i 
the potential energy for a one-component jellium model converts into 

n^'' = |tE^(^^^/^o) + ^ E '^K)|^equai(nfc)p + ^e , (US) 

with a natural extension to the general case 

n^'^'^ = li E^^^^-^(^^V^o) + ^ E '^K)I4(^.0P + ^f • (114) 

B. 2D case 

The generalization of the formulae found in a square 2D geometry to a rectangular simu- 
lation box comes in a similar manner. It is sufficient to take the resulting expressions for the 
two-dimensional problem 06 3 p and (16 7p , and to perform the necessary substitutions n ^ n^ 
and n — )■ n^, 

^^-^^ -2^ r(fc/2)|„„ + .|'= + 2-,__ fjiT^j h {-^^} + '''«.. (115) 
«-2^^ r(*/2)K|' +__Lf(I7^^|(^^J+'/v -F(|TT)' '"«' 

where V'ST" i^ given by the expressions fITT]) . f l71|) or fl75]) . 

For a long-range interaction within the jellium model, the potential energy becomes 

n^^' = 7iE^(^^./^o) + |^e (117) 



T k / ^ r \ I.J I u/ ■ r)T k 



with the notation 



Lo = {L^Lyf^ (118) 

n, = n^L^/Lo + UyLy/Lo (119) 

nfc = n^Lo/L^ + HyLo/Ly . (120) 
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For a inulticoinponent gas (quasi-neutral in case of a long-range potential), the potential 
energy is 

n'^"^ = TF E mMn,/L,) + ^e ■ (121) 



^0 i<j- ^^0 



Finally, the usual modification to calculate the momentum space sum linearly in N is 
given by 



n^'' - 7i E^(^^^/^o) + ^ E <rik)\~S.^Unk)? + ^e (122) 

n^''^= TfcE^^^^-^("^^/^o) + ^ E '^K)l4K)l' + ^e , (123) 

with ■?/', ^, S'equai, 'S'g defined by (11091) - (11121) in their corresponding two-dimensional variants. 



VI. EQUATION SUMMARY 

In the previous sections, we have derived general expressions of the Ewald sums for 
polytropic l/jrl'^ potentials in three- two- and one-dimensional systems. For integer values 
of fc, the polytropic potential reduces to a power-law interaction, which comprises realizations 
of high physical relevance. Integer power-law potentials include 

• k = 1 - Coulomb 1/|^| interaction; 

• k = 2 - Calogero-Sutherland 1/|^P interaction; 

• fc = 3 - dipole-dipole l/|r|^ interaction; 

• /c = 4, 5, 6 - interaction between different Rydberg atoms; 

• /c = 6, 12 - Van der Waals interaction. 
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The expressions for the potential energy for both the jelhum model and the general case 
of a charge-neutral simulation cell are the following 



n^''' =Tk E m^i^^j/Lo) + ^i (124) 

n^^^=7lEv^(-^^/^o) + f^e (125) 

V^(r) =^ R{nr, r) + Y, K{nu, r) + d (126) 

e =Y,(pi^'^r) + «:(nfc)) + Ci + C2 (127) 

R{n,r)=p{n + r) (128) 

/^(n, r) =/t(n) cos(27rnr) (129) 
2^4a'=-3 if A; ^3 



C?°=<J ('=-^)^[t] ^ (130) 

-An + 47r ln(2a) if A; = 3 



27ra'= 



m if M 2 



Cf =<! ('=-2)r[|] ^ (131) 

27rln(a) if A; = 2 



a 



^^ J(W,L.)V3.,3D ^^^^^ 

[ (L,Lj,)i/2 in 2D 
rir ={n ■ L)/Lo, with L = {L^, Ly, L^) (134) 

rik =(n ■ X')^o, with L' = (1/L,, 1/L„ 1/L,) . (135) 
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Alternatively, by performing a momentum space sum the above set of equations become 



ntT^O 



2U 



n 

1= Y^ p{nr) + C1 + C2 

5'cquai= g- y^exp(27rinfcrj/L) 
j 
Sg= y^ Qj exp(27rinfcrj/L) . 



(136) 

(137) 
(138) 
(139) 
(140) 
(141) 



In accordance with considerations discussed in preceding sections, the simulation cell has 
to fulfill the charge neutrality condition {Y2i=i 1% = 0) fo^ long-range potentials. Also, notice 
that in the particular case of a cubic simulation cell, 1%^ = rik = n. 

Explicit expressions of the coefficients p{n) and K,{n) for the most relevant interactions 
are summarised for 3D and 2D systems in Table 1 and Table 2, respectively. 



TABLE I. Table 1. Coefficients p(n) and K{n) taken from Eqs ^ and i^ for 3D geometry. LR 
and SR stand for long range and short range, respectively. 



p{n) 



K{n) 



LR -i— r 

\r\ 



erfc(«|ri|) 



re «^ 



LR 



w 



^erfc^^ 



SR T— rj 






27r I ^ae 



vr n errc-!— ^ 



SR -j— re 



erfc(a|ri|) , 4e~™ "' 



3V.H-^ + («H)^) 



47ra^ 



/I ft V V V V 









oJrv -j — jt; 

Vv 



(^ + fT + ;i^)^' 



SR T— TT 



E 

m=0 



(an) 



2m — a^n^ 



m! 



n 



12 



(-2r«(7-2™)!!(=^)'^"£^ 



m=0 

320F(^) . ,rn 



945 
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VII. PRACTICAL APPLICATION AND OPTIMISATIONS IN THE EWALD 
TECHNIQUE 

A. General notes 

The basic idea of the Ewald method is to calculate slowly decaying sums in a rapid manner 
by means of the Fourier transform of the slowly converging part. Although conceptually it 
provides an exact result, the number of terms which has to be summed in order to reach 
the needed convergence is a priori unknown. Once we choose the interaction potential, this 
fixes the exact form of the sums to calculate, and the practical remaining question is the 
proper choice of the free parameter a and the numbers of terms to be calculated in both 
sums: Nr and Nk in coordinate and momentum spaces, respectively. The computer time T 
is a function of only N^ and A^^, T = t^N^ + tfc^fc, with the constants t^ and tk depending 
on the complexity of the coefficients in the sums. One can notice that t^ is usually much 
less then t^, since in the Jacobi-transformed sum we only calculate cosine functions, which 
is generally far less time-consuming than the complicated functions appearing in R. It is 
clear that the parameter a affects only the resulting error in the energy. In fact, the value 
of a being very small or very large eliminates errors in one of the sums, but amplifies them 



TABLE II. Table 2. The coefficients p{n) and K{n) taken from Eqs ([37]) and ([77]) for 2D geometry. 
LR and SR stand for long range and short range, respectively. 



p(n) 



K{n) 






erfc{o|n|) 
\r^\ 






SR T— n 



2a 



oP-n^ , erfc(a|Ti|) 



+ 



H^ 



4 ( yJHae 
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in the other, so there is an "optimal" point for a, yielding a minimum error in the total 
energy. 

In the following, we discuss a way for error (SE) minimization assuming the calculation 
time T fixed. From our point of view, a useful approach for practical implementation is 
represented by the following scheme 

• We determine a time law T = trNr + tkNk in a preliminary calculation and fix the 
values of tr and tk. 

• We take a set of configurations, corresponding to the equilibrated state using an initial 
Ewald summation. Then, we calculate the exact energies E^^ (as a converged result 
of the Ewald summation) and the energies E{a, Nr, Nk) biased by a choice of Nr and 
Nk- For each pair (A^^, Nk), we find an optimal value of a = aopt(^r, ^A:)- 

• We choose the goal accuracy ^-Eacc (normally, well below the statistical error). We 
plot the error as a function of the computer time spent and choose the less time 
consumption case among the points that lie below 6Es,cc, therefore obtaining all the 
parameters required: a, N^ and Nk- From now on, these parameters are used in actual 
simulations. 

B. Example of optimisation 

Let us illustrate the scheme proposed in the preceding subsection taking as an example 
the problem of two-dimensional zero-temperature Bose gas of particles, interacting through 
the l/l^l^ potential. The model corresponds to the dipole-dipole interaction with all dipole 
moments aligned perpendicularly to the plane of motion. To describe the ground-state 
properties of the system we use the variational Monte Carlo (VMC) method and a Jastrow 
wave function with a two-body correlation factor which is solution of the two-body scattering 
problem |42|. 

The optimisation is done by averaging over A^conf = 50 uncorrelated VMC configurations, 
sampled according to the chosen probability distribution. We define the error 6E{a) as 
a sum over A^conf configurations of the difference of the Ewald energy E{iconi, Oi, N^, N^) , 
calculated for a given set of parameters (a, A''^, A^^) and the converged energy -Eex(^conf) = 
\miNk^oo^'^'^Nr^ooE{iconi,(y.,Nr,Nk). The dependence of the computer time T, needed for 
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FIG. 1. Dependence of the calculation time T on the number of terms N^. in the coordinate space 
for fixed numbers of terms in the momentum space N^ = 5, 9, 25, 45. 

the evaluation of Ewald sums, on the parameter set is shown in Figs [1] and |2l In Fig. [H 
we show the dependence of T on the number of terms Nr in real space for different fixed 
numbers of terms Nj^ in the momentum space. The computation time is proportional to the 
number of terms and the resulting dependence is linear in A^^- A fixed number of terms Nk 
requires a certain amount of calculations which results in a constant shift. Similarly, keeping 
Nr fixed and varying N^ produces a linear dependence in N^. with a constant shift which 
depends on A''^, as shown in Fig. [2l 

As one sees in Figs {!] and [2l the time dependence is linear both on Nk and N^., although 
the point corresponding to (0,0) in (A^^, N^) does not necessarily gives T = 0, since the 
reported time also contains some initializing calculations. The total error in the potential, 
as it is defined above, is given by 
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(142) 



According to our previous considerations, in the case of very small or very large values of 
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N. 



FIG. 2. Dependence of the calculation time T on the number of terms Nk in the momentum space 
for fixed numbers of terms in the coordinate space N^ = 1, 5, 9, 21. 



a the error coining from one of two sums, that is in the real or momentum space, grows 
and dominates over the error coming from the other sum; for a certain "optimal" range 
of a these two errors are of the same order. Notice that for each particular configuration, 
and each pair (A^^, Nk), it is possible to find aopt(^con/), such that E{aopt{iconf)iiconf) — 
Ecx{iconf) = 0. Instead, our task is to obtain a "universal" parameter ao, minimizing the 
total error f ll42p . The mean over the configuration set of the biased energies ^(aoj^conf) 
is used as an estimation for the mean of the exact energies E^x, introducing an inevitable 
systematic error. As it appears in typical calculations, this error is at least one order of 
magnitude smaller than the statistical error f ll42p given by the minimization of SE{a). 

A second step is the study of the dependence of the error and time on different pairs 
[Nr, Nk). The calculation time can be split as the sum of times for summing up in real and 
momentum spaces, 

T = Nrtr + Nktk (143) 

with Nr, Nk being the numbers of terms in each sum. Every one of these sums converges 
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FIG. 3. Resulting error as a function of the computer time for different parameter sets. 

when A^r, A^^ — )■ oo to a certain value, depending on a, while the sum of the limiting values 
is a constant. We can take into account the errors, corresponding to each of the sums 
separately. For a — )■ the error for the real space term is zero and the other one tends to 
infinity (and vice versa as a — > oo). The minimum total error should therefore correspond 
to the value of a, satisfying the relation d{5Er + SE^)/ da = 0. 

Focusing on the 2D system of our example, we note that the long-range expansions of 
the terms in (!37l) and ( |77l) are similar, in a sense that the leading terms in both expressions 
are Gaussians, 
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(144) 
(145) 



The power-law terms in n and the constants Cr, Ck may be neglected since the leading 
behaviour is driven by the Gaussian. The cut-off errors due to finite numbers of elements in 
the sums can be evaluated by ignoring the discrete structure of the images and approximating 
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the sums by uniform integrals, 
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(146) 



with R ~ ^jNr/'K and K ~ y^Nk/n the approximate cut-off lengths in real and momentum 
spaces, respectively. The optimal value for a can be obtained by solving the equation 
dSEr/ da = —d6Er/ da. The first-order approximation of this equation is found by taking 
logarithms of both sides and omitting constants and terms, depending on a logarithmically, 
that is 

Ak/a^ - Ar-a^ = (147) 

with Ak = T^Nk and A^. = Nr/ir, which yields 



a = {Ak/ArY^^ = (n^Nk/Nr) 



1/4 



(148) 



Then, at lowest order one finds (11461) . 



5E ~ exp{-a'^Nr/TT) = exp{-^/NkNr 



(149) 



Since the calculation time is linear with the numbers of elements Nr and A^^, we may conclude 
that with Nk fixed and comparatively large N^, \n(5E) ~ y/N^- ~ vT and vice versa, with 
Nr fixed and large N^, \ii{SE) ~ ^/N^ ~ vT- This power law may be easily checked in our 
calculations, as it is shown in Fig. ([3]). Note that for the obtained value of a the errors of the 
real- and momentum-space cutoffs are of the same order of magnitude, that is 6Er ^ SE^, 
which may serve as a rough criterium to optimise the parameter a. 

A more advanced procedure for optimisation of the parameters, proposed by Perram et 



al. 



B|, yields an asymptotic scaling N^/'^ ^ with A^ the number of particles. It is based on the 



form of Ewald summation with the momentum space sum, linear in A^ (11361) . Suppose the 
values of the calculation time tr, tk to perform unit computations in both sums are known 
and the target error level exp(— p) is fixed. Then, the total execution time in the real and 
momentum spaces is 



T = Tr + Tk = N^ttRX + N-KKHk 



(150) 



34 



with p = a^R^ = Tc'^K'^/a'^. Expressing K as K = p/{ttR) we can see that the minimum of 
the total time T corresponds to 

R,„,=|l!V'Y^V"Ar-./. (161) 



opt \ ; \ J. 

VtT/ \tr. 



--Kf)"(r)"-"' (-' 



"opt =Vvr 



-1/4 



ArV4 _ ^^53^ 



The computation time is equally divided between the real and momentum space parts (this 
was also stated in our simple optimisation scheme), with a scaling of the whole summation 
given by 

T = 2N^nRHr = 2py/WkN^^^ (154) 

Notice that the values of the free parameters change very slowly when the simulation cell is 
enlarged, and in particular a is not affected by the choice of the precision. Similar formulae 
for the optimised parameters in three-dimensional systems, with a discussion of different 
techniques to improve performance of the Ewald summation, are given by Fincham [43|. A 
more precise and detailed analytic study of the cut-off errors with verifications of the analytic 
results in actual calculations can be found in the work of Kolafa and Perram [4^ . An opti- 
mised method for treating the truncation error in Ewald sums with generic potentials was 



proposed by Natoli and Ceperley 45|. While the needed CPU time scales as 0{N\\iN)^''^ ^ 
it was shown that in the example of the Coulomb potential the method resulted in greatly 
improved accuracy compared to that of standard Ewald technique for a comparable com- 
putational effort. This method is based on an expansion of the real space function in an 
arbitrary radial basis with a parametric set of numbers in place of the fc-dependent pref- 
actors of exp(27rinr). The subsequent minimization of x^ with respect to the whole set of 
parameters yields a final optimal solution, that is the real space expansion coefficients and 
the fc-space factors. This technique was also applied to derive the optimised summation 
formulae for the two-dimensional Coulomb system 46|. 

In general, the unit computation time in momentum space is 2-4 times faster than the 
one in real space. Taking the following reasonable assumptions p = An, tk/tr = 3, we find 
i?opt ~ 2.6/N^^^. We want R to be below 0.5, since in this case the summation in the real 
space reduces to the accumulation of the single component n = 0. This condition i?opt = 0.5, 
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with our previous assumptions, corresponds to 

iV„pt = 770, iropt = 8.0, aopt = 7.1 (155) 

In smaller systems, the other components of the real sum, starting from \n\ = 1, should be 
considered. 

It is worth pointing out that if the interaction is very strong at short distances (as for the 
Lennard- Jones potential), then in principle the real-space cut-off R can be chosen below the 
"hard core radius" with a large enough value of a. This leads to the possibility of dropping 
completely the real-space part of the total sum and treat the fc-space only. This can be 
advantageous in different aspects, especially with the current progress in the development 
of efficient FFT-based methods. Nonetheless, we are not aware of any present application 
of a similar technique. 

VIII. CONCLUSIONS 

In the present work, we have applied the Ewald summation method to l/lrl'^ polytropic 
potentials in three-, two- and one- dimensional geometries in a simulation box with periodic 
boundary conditions. We have found the explicit functional forms for all the components 
of the sums in both real and momentum spaces, with special attention being paid to the 
cases of long-range interactions, that is conditionally convergent or divergent potentials 
(corresponding to k < D, with D standing for the dimensionality), "marginal" interactions 
{k = D), and short-range interactions (with k > D). For the latter case of short-range 
interaction potentials, where in principle a straightforward summation of the initial sum 
([1]) is possible, the Ewald method is shown to be useful, as it yields the faster (Gaussian) 
convergence rate. A condition of charge neutrality of the simulation cell is stated to be 
necessary for conditionally convergent and divergent potentials; a homogeneous positive 
charge background ("jellium" model) is introduced as the most relevant and frequently used 
kind of neutralization. The conditionality of the convergence for a charge-neutral system, 
governed by the Coulomb interaction, is discussed with a justification of the use of a specific 
periodicity-preserving convergence factor. The derivation technique, presented in our work, 
is consistent with the arguments of de Leeuw et al. |2|. 

The results are first presented for the case of a 3D system in a cubic simulation box in order 
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to explain the general mathematical procedure, which for the specific case of the Coulomb 
potential recovers well-known results [47|. Later on, the same mathematical technique is 
applied to 2D and ID geometries. For the one-dimensional case the initial sum for the 
potential energy is explicitly evaluated ( PB]) . nonetheless the Ewald summation is developed 
for this case too and may be used as a mathematical equality. The special representations 
of the reciprocal space sums, linear in the number of particles A^ and hence more efficient in 
actual modeling, are presented for 3D and 2D systems. The explicit expressions for the terms 
of the Ewald sums are given in a tabular form for physically relevant potentials with small 
integer power indexes k, as dipole-dipole interaction potential, Lennard- Jones potential and 
others in both three- and two-dimensional geometries (see Tables 1 and 2). 

When the simulation box cannot be chosen cubic, for example in a modeling of a three- 
dimensional hep crystal structure, the Ewald method can also be applicable after a certain 
modifications. Formally, it consists in the choice of an appropriate rectangular simulation 
box and a substitution of the vector n by n,. = {n^Lx + TiyLy + rizLz)/ Lq and n^ = 
ijix/ Lx + ny/Ly + nz/Lz)Lo in the real and momentum space sums, respectively [see fll02p 
and (IT03D . 

The optimisation of the involved parameters, that is cut-off numbers in both sums and 
the integration parameter a, is a necessary operation in order to improve the convergence 
rates and avoid excessive calculations. The main idea of the optimisation, proposed in the 
present work, is to perform a benchmark calculation, minimizing the variance of the result. 
A particular example of the application of the technique is presented for a calculation of the 
potential energy of a two-dimensional gas of dipoles, aligned perpendicular to the plane of 
motion. This practical optimisation technique is thought to be efficient for stationary and 
nearly uniform systems that appear, for instance, in Monte Carlo simulations. In spite of 
being very simple, it allows to find rather quickly adequate parameter ranges. The analytical 
estimations of the parameters are given as well and are proven to be consistent with the 
results, obtained in our method. A more sophisticated method to optimise the calculation 
parameters, taking advantage of the 0{N) representation of the Fourier transform sum, is 
also presented with explicit estimations of the parameters for a typical system simulated by 
Quantum Monte Carlo methods. 
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APPENDIX 

We prove that the sums 5*-+ and 5*++ P6] - [47|) vanish on average, allowing to calculate 
the potential energy over the negatively charged particles' positions only. 

• First, let us show that the integral of ip over the cell is zero. Since the distances are 
in the units of L, consider the cubic cell Vt = (x, y, z) G [—1/2, 1/2]^, that yields 

/ i,{r) dr = Ji + Ja + Ci (156) 

where 



It can be easily seen, that the second integral J2 is zero 



J^= I drVi?(n,r) (157) 

J^= I drVis:(n,r) . (158) 
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As far as the integral Ji is concerned, we can notice that the regions Q'{n) = r + n, 
where r ^ Q, n ^ 1? are the same cubic unit cells, displaced by an integer vector, 
thus covering all the coordinate space with only zero-measure intersections. It means 
that the summation of the integrals in (11571) over the cell ^2 can be substituted by the 
integration over the whole coordinate space, 

r(|,aV + np) 



Ji=Y.\ Riri,r)dr = J2 [ 
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and thus the whole integral (I156P is equal to zero 
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Consider two species of the particles: negative charges qi on positions r^ and a pos- 
itively charged and uniformly distributed background with a total charge q+N+ = 
—QiNi, ensuring the neutrality of the cell. Let us demonstrate that S'_+ is equal to 
zero, when the number of background charges tends to infinity. In this case the sum 
fl46p for S^^ may be rewritten as an integral over the background charges' positions 

S^.^j:^.[iir,-r.)..K-j:..fi'ir).,r, (161) 

where we did the change of variables r = Vp — r.i. The regions Vt and Vti refer to 
the original simulation cell and the same cell, moved by the vector rj, and a stands 
for the background charge density a = —qiNi/V{Q). It is clear that every vector 
r = {x, y, z) G Vti can be displaced into the cell Vt by the corresponding shift r = 
{x, y, z) = {x — aL, y — bL, z — cL) G Q with integers a, b, c. The Jacobian J of the 
change of variables r — )■ f is obviously 1. On the other hand, due to the periodicity 
of?/', 

tP(r)=tp{f), (162) 

and r runs over the whole region Q due to the conservation of the volume with J = 1. 
Finally, Eq. (I16ip can be written as 

^_+ = y'g, / ^(f)(Tdf = . (163) 

. Jn 

In the similar manner, the interaction between the charges of the background S-^-^ in 
the limit A^+ — t- oo is given by the double integral 

dr2ij{ri - r2)cr^ 
n Jn 



S++ = - I dn / dr2^(ri - r2)a' = , (164) 



since J^ip^ri — r2) dr2 = 0, following the same arguments as for the case of 5. 
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